Long-term stability of collisionless equilibrium configurations; 

algorithms 



Jean-Rene Gauthier 1,2 and Hugo Martel 2,3 



in 
o 
o 

(N 
C 

a 
>—> 

oo 

(N 



> 

o 
in 
o 

Or 

6 



ABSTRACT 

We present the summary of the theoretical aspects and algorithms used in an undergraduate 
(JRG) summer project based on numerical N-body simulations of collisionless systems. First, 
we review the importance of numerical N-body simulations in astrophysics. We introduce the 
different codes used and their performances. We then introduce four famous density profiles : 
Hernquist, NFW, truncated isothermal sphere, and Plummer. The history of these profiles and 
their dynamical properties are discussed in the third section. In the fourth section, we present 
the Barnes & Hut tree code, its features and performances. We then explain how to build and 
incorporate a multiple time stepping scheme in a tree code. The fifth section is dedicated to 
the different physical measurements we used to characterize the dynamics of the N-body system. 
Finally, we describe the future work that will be done with these codes, mainly the study of the 
adiabatic growth of a black hole at the center of a spherical distribution of stars. 

Subject headings: methods: n-body simulations — Galaxies: Kinematics and Dynamics — Galaxies: 
Structure — cosmology: Miscellaneous 



1. Introduction : numerical N-body simu- 
lations in astrophysics 

Numerical simulations, and particularly N- 
body simulations, play an important role in as- 
trophysics and cosmology. Before the advances 
in computer developments, it was impossible for 
astrophysicists to compute the evolution of a large 
system of particles evolving under their gravita- 
tional influence. The amount of calculation grow- 
ing like 0(N 2 ), the direct summation for more 
than a modest number of particles is unimagin- 
able. In 1941, Holmberg (Holmberg 1941) did a 
pioneering study of the tidal disturbances due to 
a close encounter between two nebulae. He stud- 
ied the tidal deformations by reconstructing the 
orbits of each mass elements. He used light bulbs 
to represent these elements (replacing gravitation 
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by light). The power of the bulbs is proportional 
to the mass and the total light was measured by 
a photocell. Each nebulae was represented by 37 
light bulbs. Even with this primitive but quite 
ingenious setup, the creation of spiral arms was 
detected and a good estimate of the energy loss 
was also provided. 

During the 60's and 70's with the development 
of computers, scientists developed direct summa- 
tion codes to evaluate the potential exerted on 
a few hundred particles. In 1970, Peebles stud- 
ied the collapse of a cloud of 300 particles as a 
model of cluster formation (see Klypin 1996). Six 
years later, White studied the collapse of 700 par- 
ticles with different masses. Unfortunately, with 
these low resolution simulations, numerical arti- 
facts (two-body relaxation) are important and the 
results from these simulations were not reliable. 
The first cosmological simulations were made in 
the mid 70 's. Most of these simulations of struc- 
ture formation confirm the existence of a flat spec- 
trum of initial fluctuations. 

The last thirty years have seen the development 
of several N-body codes. Each one is particularly 
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suitable to study a specific range of dynamical 
problems. Each algorithm has its strengths and 
weaknesses. In the next few paragraphs, we de- 
scribe the characteristics of a few popular codes 
for collisionlcss systems, i.e. the study of internal 
dynamics of galaxy, interactions between galaxies 
and the clustering in an expanding universe (Sell- 
wood 1987). 

1.1. Direct summation codes or particle- 
particle method 

In this algorithm, the force is computed ex- 
actly for each particle. If the system consists of 
N particles, the number of operations to evaluate 
the force exerted on these particles is proportional 
to N 2 . This method is flexible, but has a high 
computational cost and is limited to systems with 
N rts 1000. Usually, the system is integrated us- 
ing a second-order scheme like leapfrog or Runge- 
Kutta. The developers of these codes were aware 
of the fact that when particles were closed to each 
other, the acceleration becomes important and one 
might use multiple timesteps (see section 3) in the 
integration scheme in order to compute the parti- 
cles trajectories accurately. 

1.2. Grid Methods (Particle-Mesh) 

This method is the fastest when one has to deal 
with large N. The number of computation is of 
order 0(N + N g \ogN g ) where N g represents the 
number of grid points. This method determines 
the force acting on each particle by evaluating the 
gravitational field. The trick is to associate parti- 
cles to nearby point mesh. Hence, each mesh point 
has a defined density attribute. The gravitational 
potential is evaluated by solving the Poisson equa- 
tion 

V 2 = 4nGp (1) 

on a number of grid points within a fixed volume 
in space. The force on each particle is obtained 
via interpolation techniques. 

There are several ways to assigning particles to 
the mesh, but one must be careful with the fluctu- 
ations of the force when the particles are close to 
each other. The continuity of the derivatives the 
function used to assign particles to mesh is the 
first criterion. 

The simplest assignment scheme is the "Nearest- 
Grid-Point" (NGP) scheme. The density at a 



mesh point is determined by the number of par- 
ticles in the cell centered on the point divided 
by the volume of the cell. However, the forces 
are discontinuous and this technique is rarely 
used. There are other more accurate assignment 
schemes. Two of them are the "Cloud-in-Cell" 
(CIC) and "Triangular-Shaped-Cloud" (TSC) . 
The former one gives a continuous force, but the 
first derivative is discontinuous. The latter is more 
accurate. It employs an interpolation function 
that is pieccwise quadratic and each particle is 
assigned to a larger number of mesh points. 

The grid methods have been successful in sim- 
ulations of clustering in an expanding universe 
and in studying galactic disk dynamics. However, 
these techniques have several limitations. First, 
the spatial resolution is roughly the distance be- 
tween mesh points. The geometry imposed by the 
grid cannot be adjusted to the changing shape of 
the system. Also, the evaluation of the potential 
on grid points where there are no nearby parti- 
cles is a waste of CPU time and a weakness of the 
method. 

1.3. Tree Codes 

The goal of a tree code is to regroup parti- 
cles to reduce the number of contribution to the 
force acting on any particle. Basically, this tech- 
nique has a lower computational cost than the di- 
rect summation scheme. The calculations grow 
like 0(N log N). In this scheme, the contribu- 
tions from nearby particles have to be summed 
directly, but the particles that are far away can 
be regrouped and replaced by a monopole term in 
the summation. At each step, the "grouping" can 
be different because the particles move in space. 
Hence, one must rebuild the tree at each time step. 
This is one drawback of the method. The first two 
well-known tree codes were designed by Appel in 
1985 and Barnes & Hut in 1986. Section 3 is ded- 
icated to the Barnes & Hut tree code. 

This technique has advantages and weaknesses. 
First, the tree structure takes into account the 
shape of the system (Barnes & Hut 1986). There is 
no fixed geometry. Second, there is no time lost in 
evaluating the potential in regions devoid of mat- 
ter. However, this technique is slower than PM 
methods. Also, it is more difficult to implement. 

Finally, the method has been particularly suc- 
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cessful in studying the collisions between galaxies, 
the internal galactic structures and cosmological 
structures formation. 

2. Models 

The previous paragraphs introduced basic con- 
cepts related to numerical iV-body techniques. In 
this section, we present the main physical quanti- 
ties of four different spherical models used and/or 
discovered in numerical simulations : Hernquist, 
NFW, Plummer, and TIS density profiles. For 
each one of them, we include the cumulative mass 
function, the gravitational potential and the distri- 
bution function if an approximation or analytical 
expression can be found. 

2.1. The Hernquist density profile 

Hernquist (1990) first introduces this density 
profile as a fit to the de Vaucouleurs surface bright- 
ness profile, 
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where I(R e ) and R e represent respectively half of 
the total luminosity and the isophote value cor- 
responding to it. The de Vaucouleurs profile is a 
good fit for the surface brightness of most elliptical 
galaxies. It is also a good approximation for the 
distribution of stars in the bulge of spiral galaxies. 
The profile has the following shape, 
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M is the total mass of the system and a is a scale 
length. The potential can be found via the Poisson 
equation. Using the preceding equation, one finds 

<f>(r) = — . (4) 

r + a 

It is possible to express the density in function of 
the potential <j>(p) This enables us to evaluate an- 
alytically the distribution function (DF) (Binney 
& Tremaine 1987) of the system, 
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The solution for the Hernquist 's profile is 
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Useful approximations can be made on this distri- 
bution to find the limiting case when q — > and 
q — ► 1. From equation (3) and (4) one can find the 
the velocity dispersion by solving the Jeans equa- 
tion. For an isotropic system [oy = ag = a(<p)] 
where r -C a, 
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The model provides a good fit for the de Vau- 
couleurs i? 1 ' 4 law [see Fig. 4 in Hernquist (1990)]. 
Its distribution function [eq. (6)] is found analyti- 
cally. 

2.2. The Plummer profile : a model for 
the distribution of stars in globular 
clusters 

The Plummer model (Plummer 1911) has been 
widely used in numerical simulations of star clus- 
ter dynamics (e.g. Aarseth, Henon & Wielen 
1974). It is a polytrope of index 5 given by : 

3 M 1 

' (r) = 4^[l + (r/*)^- (10) 

where M represents the total mass of the cluster 
and R is a scale length. The mass enclosed in a 
radius r is determined by the following equation : 

M(r) - M (11) 

i? 3 [1 + (r/R) 2 f 2 

and the gravitational potential is given by 

GM 1 
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The density profile can be expressed in terms of 
the gravitational potential. The distribution func- 
tion [see eq. (5)] can be found analytically : 
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The Hcrnquist and Plummcr models can be eas- 
ily used as initial conditions for N-body simula- 
tions. Aarseth, Henon & Wielen (1974) used a 
simple method when it is possible to express ana- 
lytically the radius in function of the mass : 



M(r) -» r(M) 



(14) 



The technique needs the creation of a subrou- 
tine which generates five random numbers X\, X2, 
A3, X4 and X5 (uniform deviate between and 
1). Assume that M(r — ► 00) is normalized and 
G = M = R = 1, then rj (the radius of the j th 
particle) is determined by rj(Xi). In the case of 
the Plummcr model, 
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The Xj, yj and zj components can be evaluated 
using X 2 and X 3 using the following equations : 



j> 3 
Vi 



= (l-2X 2 )rj 

= (r 2 - z 2 ) 1 ' 2 cos(27rA 3 ) 



(16) 
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For the velocity, one has to evaluate the escape 
velocity (V e ) at each rj. We can write Vj/V e = 
qj and insert this ratio into equation (13). The 
probability distribution of qj is then proportional 
to 

g( qj ) = q 2 (l-q 2 )y 2 (17) 

(note : the factor q 2 comes from the fact that we 
use the speed distribution instead of the velocity 
component distribution). The values of qj range 
between and 1 and g(qj) is always less than 0.1. 
By using X4 and X 5 , we can set qj = X4 if X4 and 
X 5 fulfill simultaneously the following criterion : 



0.1X 5 < g(Xi) 



(18) 



The velocity components are determined using the 
same trick as for the position. 

2.3. NFW profile as a result of cosmolog- 
ical simulations 

The Navarro-Frcnk-Whitc (hereafter NFW) 
profile (Navarro, Frcnk & White 1995, 1996, 1997) 
is the profile resulting from the formation of cold 
dark matter halos in a hierarchical clustering uni- 
verse. NFW find that the halo profiles have the 



same shapes, independent of the halo mass, initial 
density fluctuation spectrum and the values of the 
cosmological parameters. 

The profile has the following shape (Navarro, 
Frenk & White 1997) : 



P(r) = $c 

Pcrit r/r s (1 + r/r s f 



(19) 



where r s is a scale radius, 5 C is a dimcnsionless 
density and p cr n = 3H 2 /8irG corresponding to 
the critical density for a closed universe. 

The value of S c is strongly correlated with the 
value of the halo mass. Less massive halos have 
higher S c , indicating a higher rcdshift of collapse. 

From the value of the halo "concentration", 
c = r 2 ooAs (^200 is the halo radius enclosing a 
density equals to 200p cr i t ), it is possible to estab- 
lish a relation between c and S c : 

(20) 
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The velocity curve is given by : 

V c (r) _ 1 ln(l + ex) - cx/(l + cx) 
V200 ~ x ln(l + c)-c/(l + c) 

and the cumulative mass function is 
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The gravitational potential can be found using 
equation (2-22) in Binney & Tremaine and is equal 
to 

GM(r) 4irG6 cPcrit r 3 s 



$(r) 



(23) 



where M(r) is given by equation (22). 

The next step is to evaluate the distribution 
function of the NFW profile using equation (5). 
The calculation of dp/d<& shows that p cannot be 
expressed as a function of <&, the gravitational po- 
tential : 



dp _ 8 c p crit r 3 (3r 
rf$ ~ GM{r)(r + r s ) 3 



r.) 



(24) 



For NFW, equation (5) cannot be solved ana- 
lytically, we need to use numerical methods. How- 
ever, some authors (e.g. Widrow 2000) devel- 
oped useful approximations for f(E). Widrow 
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(2000) proposed the following fitting formula for 
the isotropic DF : 



F(e) = F e 3 / 2 (1 - e)" 



1-e 



(25) 



where e = —(E - Q^) / 1 AitG5 c p a - lt r 2 s , A, F , and 
q are the fitting formula parameters for isotropic 
models. P is a polynomial introduced to improve 
the fit (P = S iPj e 4 ) In the case of NFW, F = 
9.1968 x 10~ 2 , q = -2.7419, A = 5/2, Pl = 0.3620, 
p 2 = -0.5639, p 3 = -0.0859, and p A = -0.4912. 

In Figure 1, we give a comparison between the 
Hernquist, Plummer, and NFW profile for the 
mass function, density profile, and gravitational 
potential. 

2.4. The Truncated Isothermal Sphere 
(TIS) profile 

The TIS (Truncated Isothermal Sphere) is a so- 
lution of the Lane-Emden equation 



d 2 uj 2 duj 



+ T-rr = e 



(26) 



where u> = ln(p/p n ), £ = r/r and r is the core 
radius : 

•» - ■ (27) 

The TIS model corresponds to the outcome of 
the collapse and virialization of a top-hat den- 
sity perturbation. This model involves a non- 
singular, truncated isothermal sphere (Shapiro & 
Iliev 1999). For a given value of boundary pres- 
sure pt and mass Mo, there is a unique value of 
£t (where £ t = r t /ro - r t is the truncation ra- 
dius) which minimizes the total energy. This min- 
imum energy solution is the unique TIS solution 
preferred in nature as the outcome and virializa- 
tion of a sphere in the presence of a fixed external 
pressure. The virialized object has the same total 
energy as the top-hat before the collapse. Shapiro 
& Iliev found that this solution implies £ = 29.4. 

There is a useful fitting formula to equation (26) 
that is a good approximation within several core- 
radii (Natarajan & Lynden-Bell 1997) : 



p(0 = 



B 



(28) 



a 2 +e p+e ' 

For the particular TIS case, 
(A, a 2 , B, b 2 ) = (21.38, 9.08, 19.81, 14.62) . (29) 



The TIS profile seems to be a better alternative 
to the NFW profile in many cases. Shapiro & Iliev 
(2000) showed that the projected mass density of 
the cluster CI 0024+1654, determined by strong 
gravitational lensing, is well fitted by a TIS sphere. 
They also found that a central cuspy NFW fitting 
profile implies a velocity dispersion which is too 
large by a factor of <~ 2 to be consistent with the 
measured velocity dispersion. It also appears that 
the TIS model can be a good fit to the density pro- 
file of dark matter-dominated dwarf galaxies (Iliev 
k Shapiro 2001). Finally, the TIS solution for viri- 
alized cosmological halos reproduces fairly well the 
po —<J V correlation (Shapiro & Iliev 2002). In fact, 
current data suggest that the central mass densi- 
ties po of cosmological halos in the universe are 
correlated with their velocity dispersion a v over 
a very wide range of a v (Shapiro & Iliev 2002). 
For all these reasons, it is important to study this 
density profile. 

The cumulative mass function for equation (28) 
(Natarajan & Lynden-Bell 1997) is 
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and the gravitational potential is 
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Equation (28) can be integrated to yield an an- 
alytical fitting formula for the TIS rotation curve 
(Iliev & Shapiro 2001) given by 



v(r) 
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(32) 
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In summary, we have briefly introduced few ba- 
sic models used in current numerical simulations of 
structure formation, galactic dynamics, and star 
clusters. In the next section, we review the BH 
tree code, the TV-body code we used to make these 
systems evolve. 

3. Barnes &: Hut Tree code 

The Barnes and Hut (hereafter BH) tree code 
(Barnes & Hut 1986) is the one we use for the 
simulations. However, we modified its structure 
to include a multiple timesteps integration scheme 
(see section 4). In this section we describe how 
the tree code works : its main features and basic 
related concepts. 

3.1. The tree building 

The BH technique is based on a hierarchical 
division of space into cubic cells. At the beginning 
of the simulation, the particles are distributed, one 
by one, into the root cell. If two particles are in 
the same cell, then this cell is divided into eight 
"daughter" cells (an oct tree) of equal volume (see 
Figs. 1 and 2 in Barnes & Hut 1986). The aim 
of this process is to make sure that at the end of 
tree building step, two particles cannot reside in 
the same cell. 

For each cell with subcells in it, there is an as- 
sociated pseudo-particle which contains the total 
mass in the cell and located at the center-of-mass 
of the particles distribution in the cell. 

This building step must be repeated at each 
timestep because the particles move and they are 
no longer in the same cell. It takes usually less 
than 10 % of the overall CPU time to build the 
tree (Hernquist 1987). 

3.2. The force evaluation 

The significant difference with direct summa- 
tion methods (see section 1.1) resides in the ap- 
proximation made in the force calculation. We 
first consider the interaction between a group of 
particles (represented by a cell of higher level) and 
a particle. We want to use a criterion based on 
the size of the cell (/) and the distance between 
the center of mass of this cell and the particle (d) 
to verify if we can consider this group as a point- 
mass particle or if we need to go down to a lower 



level of the hierarchy and resolve its constituents. 
This criterion can be written as 



9 > 



(33) 



where 9 is a free parameter called the "opening 
angle" parameter. Therefore, to compute the force 
on particle p, one needs to walk down the tree, 
starting by the bigger cells and check for each cell 
if equation (33) is verified. If it is the case, then 
add the contribution of the pseudo-particle of this 
cell to the force exerted on p. The result of this 
"walking" process is to reduce the number of terms 
in the force calculation on one particle from N to 
log TV. 

The value of 9 is chosen at the beginning of the 
simulation. Hernquist (1987) mentioned that for 
a large 6 (9 > 1.3), the BH code can unintention- 
ally include particle self-acceleration by not forc- 
ing subdivision of cells containing the particle. It 
is worth noting that decreasing 9 from 0.9 to 0.6 
for 10, 000 particles corresponds to increasing the 
CPU time by a factor <~ 2 — 2.5. The range of 
typical values for 9 is between 0.7 and 1.0. This 
is a good compromise between accuracy and per- 
formance. The typical error for 9 = 0.9 relative to 
a direct sum is ~ 1% (Hernquist 1987). The error 
for a fixed 9 decreases as N increases. 

After establishing the opening angle criterion 
and determining which cells need to be resolved, 
one needs to calculate the force exerted by either 
cells or close particles. The particles are repre- 
sented by point masses. However, the force that 
two point masses exerts on each other becomes 
large at small distances. Therefore, the velocities 
will change very rapidly and very short timesteps 
are needed in order to properly determine the mo- 
tion of the particles (see section 4) . We can avoid 
this problem if the particles are extended (Binney 
& Tremaine 1987). To compute the force exerted 
by particle (or pseudo-particle) j on particle i, one 
can use the following formula : 



Grriimj (xj — x*i) 



(34) 



where e is the softening length. The maximum 
force occurs at \xi — Xj\ 2 = e 2 /2. This potential 
is called the softened point-mass potential (Binney 
& Tremaine 1987) also referred as plummcr-type 
softening. 
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Finally, for 8 — 1.0 and a Plummcr model (sec 
section 2.3), the tree walk process and force sum- 
mation take respectively 80 % and 10 % of the 
total CPU time (Hcrnquist 1987). 

3.3. The integration scheme 

Once the force has been evaluated for each par- 
ticle, we need to update the position and velocity 
of the particles. We describe here two different 
schemes of order 2 : leapfrog, and Runge-Kutta. 
The leapfrog technique is used in the BH tree code. 

3.3.1. The leapfrog integrator 

Suppose that we know the particles positions 
and velocities at step n and we want to advance 
the system to step n + 1. This can be done using 
a time-centered leapfrog : 

- At 

v i,n+l/2 — v i,n + a i,n~^~i 

n,n+i = n,„ + v i>n+1 / 2 At, (35) 
_ At 

Vi,n+1 — v i,n+l/2 + a n+l~^- ■ 

At is the timestep. The force is computed at the 
beginning of the timestep and is used to extrapo- 
late the position and velocity at respectively n + 1 
and n + 1/2. 

3.3.2. The Runge-Kutta method of order 2 

This method is similar to the Leapfrog scheme. 
However, it requires the storage of an additional 
variable, namely fi tn+ i/ 2 . The force is determined 
at the center of the timestep. The position and 
velocity at step n + 1 are extrapolated from the 
value of the force at the middle of the timestep 
(Fortin 2001) : 



Vi 



- At 

,n+l/2 — v i,n + a i,n~^-, 

- At 



r i,n+l/2 ~ r i,n + v i,n~^- , 
Vi,n+1 = V itn + a itn+ i/ 2 At, 
r»,n+l = n,n +V lin+ i/ 2 At . 



(36) 



Both methods give similar results. 

The BH tree code provides a technique for han- 
dling a large number of long-range interactions 
and concentrating the CPU time allowed on lo- 
cal interactions where more precision is needed. In 



this iterative process, the number of mathematical 
operations grows as NlogN. 

4. A multi timesteps scheme 

In the previous section, we described the BH 
tree code features. Now, we introduce a multi 
timesteps scheme which can be inserted in a tree 
code algorithm. This scheme is especially useful 
for studies of dynamical properties of halos con- 
taining black holes, SPH simulations of star for- 
mation, etc. The basic idea is that far away parti- 
cles evolve less rapidly than particles closer to the 
dynamical center of the system. The motion of 
distant particles can be integrated with a longer 
timestep. 

4.1. A measure of the integration time for 
each particle 

The first thing we need to do is to evaluate the 
integration time for each particle. One of the most 
commonly used criterion is to associate the veloc- 
ity and acceleration of the particle to its softening 
length : 



Ti = [3 min 



1/2' 



(37) 



where rj, ej, Vi and ttj are respectively the integra- 
tion time, softening length, velocity and acceler- 
ation of particle i. (3 is a free parameter varying 
between and 1. This is the most straightforward 
method to estimate the integration time needed 
for each particle. However, we are looking for a 
more meaningful criterion based on the orbital pe- 
riod of the particles. 

Equation (37) provides a value of r for each 
particle. The next step is to sort those value of n 
(i e 1, . . . , N) and find T max and T m j n . We then 
distribute these particles into different bins, hav- 
ing different integration timestep. The steps are 
separated by a factor 2 in time. If n represents the 
number of bin, then : 



n = 0.69314718 log 



(38) 



(n is round to the next integer). We have now a 
set of n different timesteps for the ./V particles. We 
will use the following notation for the timesteps : 



Atj, where j G 0, . . . , n — 1 



(39) 
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where Atj represents the timestep of bin number 
j. The longest timestep, namely Aibasic corre- 
sponds to Aio and the shortest one (Atf un d) is 
At„_i. Figure 2 illustrates the distribution of 
timesteps for n = 4. 

4.2. Evaluating the force acting on each 
particle 

The goal of introducing a multi timesteps 
scheme is to reduce the number of force evalua- 
tions for the particles having a longer timestep. To 
integrate the motion of the particles, the Barnes & 
Hut tree code uses a leapfrog integrator in which 
the forces are evaluated at each half-timestep [see 
previous section, eq. (35)]. In this multi timesteps 
scheme, the timesteps are variable between par- 
ticles and therefore, the evaluation of the force 
on the particles is not necessarily synchronized 
with the evaluation for the other ones. In Fig- 
ure 2, we show by red dots, when the force must 
be evaluated in a system of n = 4 bins. At the 
beginning of the simulation (t = 0), all the parti- 
cles are synchronized. Then, the code updates the 
position and velocity of each particles (using the 
leapfrog scheme) with an interval of A£f unc j/2. At 
t = Aif un( j/2, only the particles of bin number 3 
need a force evaluation and at t = 2 x Aif unc j/2, 
particles of bin 3 and 2 need a force evaluation. 
There is a simple criterion to check if the bin num- 
ber j needs a force evaluation. Let Q represents 
the number of Aif un d/2 elapses since the begin- 
ning of the simulations, i.e. t = QAtf un( j/2. The 
j th bin needs a force evaluation if 

Q mod = . (40) 

It is preferable to use integers in equation (40) 
instead of directly comparing t and Atj because of 
the possible truncation errors modifying the value 
of t during the simulation. 

4.3. Particle moving from bin j to bin i 

During the simulation, it is highly probable that 
many particles will move from one bin to another 
after evaluating their integration time. For exam- 
ple, a particle on a very eccentric orbit will need to 
diminish its timestep during the closest approach 
and increase it near the pericenter. The multi 
timesteps scheme should allow particles to move 



from one bin to another, but not at any moment. 
There is a simple criterion to verify if a particle ini- 
tially in bin j can move to a bin i. First, both bins 
(i and j) must be synchronized, i.e. their internal 
clock should indicate the same time and this time 
should correspond to the number of fundamental 
timesteps (Q/2) elapses since the beginning of the 
simulation. 

There is a simple way to evaluate the number 
of fundamental timesteps elapses for each bin of 
a set of n bins. Suppose that Pj represents the 
number of fundamental timesteps of the j th bin. 
Pj will be equal to Q if 

Q mod -^L. = . (41) 

If the following expression is false then Pj keeps its 
latest value when equation (41) was verified and 
the particle stays in the same bin. 

P = P 3 = Q . (42) 

You have to wait until equation (42) is verified 
before moving a particle from one bin to another. 

4.4. The value of n changes 

What happens if a particle of the fundamental 
bin (in our case, bin n) goes from bin n to bin 
n + 1? Following the previous argument, this can 
happen at each fundamental timestep. In fact, the 
fundamental timestep changes. It goes from At n 
to At„ + i. We must then take that into account 
and change the number of P for each bin. In this 
example, 

Pj = 2Pj for j e 0, . . . , n (43) 

and P n +i will be equal to the previous value of 
P n multiplied by 2. Q will also be modified, going 
from Q to 2Q and Atf unc j = Atf un( j/2. 

The opposite phenomenon happens if all the 
particles of bin n move to bin n — 1, i.e. we lose a 
bin. The process is exactly the same except that 
we divide by 2 instead of multiplying by 2. 

Finally, at the end of a basic timestep, i.e. if 

t = mAibasic (44) 

where m is an integer, we allow all the particles 
to change from one bin to another after evaluat- 
ing their integration time. Figure 2 (upper panel) 
illustrates the different situations examined in the 
previous paragraphs. 
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5. The measures, observations, dynamical 
properties of the system 

In the previous three sections, we describe the 
theoretical aspects of the models we want to study 
(section 2), the tree code used (section 3) and a 
multi timesteps scheme (section 4) implemented in 
a leapfrog integrator [see eq. (35)]. We are inter- 
ested in determining the evolution of the physical 
properties of the simulated systems. This section 
contains a brief description of the different prop- 
erties of the system one can measure using the 
N-body results. 

5.1. Conservation of energy 

First, the most important thing to measure 
when doing N-body experiments is the total en- 
ergy of the whole system. The energy must be con- 
served. The state-of-the-art simulations are able 
to conserve total energy within 0.1 — 0.2%. If it is 
not well conserved this is maybe due to a wrong 
choice of timesteps. The energy is calculated using 
the following formula : 



N N 



1 2 GnijiTii 



i=l j^i 



Tin 



(45) 



where is the distance between particle i and j. 

5.2. Density profile 

The evaluation of the density profile is primor- 
dial to study the dynamical properties of a system. 
In particular, if one studies the collisions between 
galaxies harboring a black hole, the profile can give 
serious insight about the orbits distribution in the 
inner few kpc. 

There are basically two ways of measuring the 
density profile for spherically symmetric systems. 
The first one is two divide the space into spher- 
ical bins of width Ar and count the number of 
particles in each bin. The bin width can increase 
logarithmically as r becomes large in order to have 
a good resolution in the center and a lower one in 
the outer parts of the system where nothing im- 
portant happens. In this picture, the density is 
evaluated using the following equation : 



Pk 



Aitr 2 k Ar k 



(46) 



where pu, ffc and Ar& are respectively the density, 
radius and width of bin k; Nj represents the num- 
ber of particles of mass nij in bin k. This method 
has drawbacks when the number of particle is not 
sufficient to resolve the inner bins. Some bins may 
be empty and others can contain not enough par- 
ticles. 

The second method consists of sorting the par- 
ticles by radius, then divide them by slices of 100 
or 1000 particles, for instance, and evaluating the 
average radius, width and density of each slice of 
particles. This method avoids low resolution bins 
but one needs to have enough particles to get slices 
thin enough for a good resolution. 

5.3. R w (t), R 50 (t) and R 90 (t) 

In the previous paragraphs, we have seen how 
to compute the density profile. We are now inter- 
ested in its evolution. Ri (t), i?5o(£) and -Rgo(i) 
represent the radii which contains 10 %, 50% and 
90% of the system mass, respectively. The mea- 
surement of these quantities gives a rough estimate 
of the evolution of the density profile. 

5.4. Free-fall time 

In the simulations, time is measured in dimen- 
sionlcss units. A good way to characterize the time 
scales is to express them in terms of the free-fall 
time 



3tt 
32Gp 



(47) 



This quantity represents the time needed for a ho- 
mogeneous sphere of pressureless material of den- 
sity p (where p can be evaluated at R = R50) re- 
leased from rest at t = to collapse to a point at 
t = t s (Binney & Trcmaine 1987). 

5.5. Quadrupole tensor evaluation 

To estimate the size of the cluster, one has to 
evaluate the components of the quadrupole tensor 
defined as 



Qij = E mfc ( rfc )* ( rfe )j 



(48) 



where mu and are the mass and position rela- 
tive to particle k. The indicies i and j account for 
the vectors components. Once the tensor compo- 
nents have been determined, one has to compute 
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the eigenvalues Q\, Q2 and Q3. The "semi-major 
axes" of the system are evaluated using the follow- 
ing formula : 



5Q t 
Mt„, 



1/2 



for i e 1,2,3 . (49) 



The ratios of the different a, give a measure of the 
system sphericity. If a\/d2 = (12/(13 = 1, then the 
system is perfectly spherical. 

5.6. Velocity dispersion, anisotropy pa- 
rameter 

The velocity dispersion for a set of particles is 
defined as follows: 



IE 



[(Vk)i 



7,12 



N 



(50) 



where Vi is the average i-component of the velocity 
over the distribution and N is the number of par- 
ticles. For a spherical system, it is useful to find 
oy, <Jg and To do this, one needs to find v r , v$, 
and v<f, These three quantities can be related to r, 
x, y, z, v x , v y , and v z using the following relations 



1 / 

- (xv x + yv y + zv z ) , 



(51) 



\Jx 2 + y 2 



— (xv x + yv y + zv z ) - v z 
(xv y -yv x ) . 



If ay = oe = (7(p then the system is isotropic. 

Another way to measure the degree of isotropy 
of a system is to evaluate the anisotropy parameter 
(3 defined as : 



(2«2) 



(52) 



where v t is the tangential velocity. If the distribu- 
tion of orbits is isotropic, then (3 = 0. If the orbits 
become radial, then (3 — > 1. For a set of perfectly 
circular orbits, (3 — > — 00. 

6. Conclusion 

In conclusion, this short review is an tentative 
introduction to the basic concepts related to nu- 
merical iV-body techniques in astrophysics. We 



introduced different codes used, physical models 
used and developed in numerical simulations, al- 
gorithms and measurements. We hope this doc- 
ument will be a useful synthesis for any beginner 
interested in numerical astrophysics. 

7. Future Work 

This document presents the various basic algo- 
rithms needed to do a study of collisionless equi- 
librium systems. Using the techniques presented 
in this paper, we will study the long-term evolu- 
tion of density profiles resulting from the mergers 
between halos containing black holes. 

We acknowledge the Natural Science and Engi- 
neering Research Council of Canada (NSERC) for 
an undergraduate summer fellowship (JRG) and 
thank Laurent Drissen and John Dubinski for use- 
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thanks the Canada Research Chair program for 
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A. Appendix material 

This appendix is dedicated to the description of an algorithm evaluating the number of possible two-body 
bound systems which can form in iV-body simulations of dense star clusters. The cluster has to be dense 
enough so that the softening length is short to allow the formation of bound subsystems. The following 
algorithm is a tentative method to identify possible two-body systems in TV-body files. Those systems can 
be referred to as "binary stars" . 

A.l. Brief Overview 

The Binary Stars Finder Code is an analysis tool that can identify the presence of binary systems in a 
dense cluster of particles. To work properly, the code needs to receive as inputs the dynamical data of the 
system, i.e. the position, velocity and acceleration of each particle. Usually this can be done by using typical 
output files from N-body codes. To be able to identify correctly the formation of potential binary systems, 
we need to set up conditions for the formation of such systems. 

The first one is related to the distance between particles. Suppose that particles A and B are constituents 
of a binary system. Our first hypothesis is that : A is the closest neighbor to B and vice-versa. Hence, if 
there is a third body in the very close vicinity of A, even closer than B, we can not consider the system 
formed by A and B to be a binary one. The probability that A will interact strongly with this intruder is 
higher than with B. The condition of reciprocal closest neighbor must be fulfilled. 

Second, a cluster consisting of several particles will remain bound if the total energy has a negative value. 
This is also true for binaries. So we have to sum the energy of the two members. If the results respects this 
second criterion the likelihood that we have identified a real binary increases. 

The last condition consists of giving a certain level of "quality" to every system who already satisfies the 
first two criteria. This can be achieved by determining the number of bodies that can potentially disturb 
the path of one or even two members of the binary systems. If these intruders are close enough to the 
considered system, they can break up the whole system in the next few iterations. In section A. 2. 3, we 
describe a method to reduce the problem to a single moving particle that follows a precised path around a 
fixed partner. The characteristics of this orbit will enable us to set up a "neighborhood of influence" and 
after that to count the number of "bad" neighbors. 

Now the question is how does it work ? This is the object of the next section. In fact, particular features 
and algorithms will be discussed. Also, in the last section of this presentation document, We suggest several 
improvements that could be made in further versions of the code. 

A. 2. Principal Features 

A. 2.1. "Divide- and- conquer" distribution method 

As we saw in the previous section, the first criterion consists of a reciprocal closest neighbors selection. To 
do this, we can simply compute the distance between particle X and every other particle of the cluster and 
find the closest neighbor to particle X. We could repeat the process for every member of the cluster. Since 
this calculation grows like N 2 , for a typical cluster of 500, 000 that represents a total number of 2.5 x 10 11 
mathematical operations. This method has a very high computational cost and should be avoided if possible. 

Instead of doing a direct calculation between each pair of particle, we can divide the volume that contains 
the whole system into individual cubic cells. After that, the comparisons could be made under the assumption 
that members of binary systems will be close to each other, i.e. members of the same or adjacent cells. Hence, 
if we use this method the number of operations (S) will grow approximately like : 

N 2 

E oc — (Al) 
n 

where N represents the number of particles in the simulation and n the number of cells. This method is 
called "divide-and-conquer." Doing it this way we can save precious CPU time. 



12 



The side length of each cell should be chosen so that it is greater than a few times the average distance 
between particles (to have the closest neighbor in the same cell or in an adjacent one). In order to have a 
significant gain in CPU time, the size of the cell should be smaller than the size of the whole cluster. 

Now what happens if the closest neighbor is not in the same cell. This could be possible if the particle 
(particle A) on which we want to identify its closest neighbor is near the walls of the cell. We can handle 
this kind of situations by comparing the distance between A and the closest neighbor of A in the same cell 
(value called d) with the distance between A and the walls of the cell (k, where i goes from 1 to 6 - a cube 
has six faces). Now if 

d>k (A2) 

we must evaluate the distance for particles in the cell labelled i. The process is repeated for adjacent cells 
who agree to condition (A2). Figure 3 illustrates the previous considerations. We can repeat the same 
process for each one of the 26 adjacent cells of the one we consider. 

A. 2. 2. Evaluating the Energy of the Bound System 

The second step in the identification of binary systems consists of computing the total energy of the 
two-body system. We can compute the total energy of the bound system by using the following formula: 

£tot = 7^1 N 2 + 7; m i\ v 2\ - i z * (A3) 
2 2 |ri-r 2 | 

where the subscripts 1 and 2 denote body number 1 and 2 respectively. Of course, the system will remain 
bound if E tot has a negative value. 

A.2.3. Number of neighbors in the close vicinity of the two-body system 

The last step in the identification of binaries consists of evaluating the number of neighbors in the vicinity 
of the system. First, we must characterized the size of the system we are studying. Figure 5 is a schematic 
representation of a typical two-body system. R represents the center-of-mass position vector. We can reduce 
this system to a single particle moving around a fixed massive particle. The motion of this particle must 
obey to the specific condition: its angular momentum and energy must be the same as for the previous 
system and should be conserved along the path. We can evaluate the reduced angular momentum using this 
formula : 

\Ltot\ = \Li\ + \L 2 \ = itrv = fiy/GMa(l-e 2 ) (A4) 

where f = f 2 — rj a is the semimajor axis of the elliptic path of the moving particle and \x is the reduced 
mass of the two-body system: 

H = ■ . A5) 

mi + 777-2 

For a single particle moving around a fixed one, the energy can be written as : 



and using the virial theorem: 



E=\^-^t (AO) 

E .^-°M>.-<%L. (A7) 

By replacing a in equation (A4) by its value in equation (A7), we can find the ellipse eccentricity e : 

p — I 2 |£tot| 2 |£tot| (M) 
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Once the eccentricity has been evaluated, we can set the "influence" radius to be the aphelion distance to 
the center-of-mass of the system. The aphelion radius is determined using the following : 



So, if the distance between a particle i and the center of mass of the system is shorter than (3r a we can 
consider this body to be in the close vicinity (C.V.) of the system j: 



where (3 is a free parameter. It will probably disturb the path of one or maybe two particles. 

A. 3. Further Improvements 

Several improvements will be made in a next version of the code. Currently, the code can only handle 
systems contained in a cubic volume. Of course, we could put the whole cluster in an augmented cubic 
volume but there will be many empty cells and this is not really optimized. Modifying the code so that 
systems with rectangular shape can be well-treated is a first thing to do. 

The current algorithm analyzes only one snapshot data of the system and evaluates the formation of 
binary systems. It could be really interesting if the code could integrate the motion of binaries with the data 
of several snapshots taken at different times. By doing this, we could be able to tell if a binary we identified 
previously is going to break up or not. 



r a = o(l + e) 



(A9) 



(A10) 
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Fig. 1. — Density, mass and gravitational potential for three models: solid line: Plummer model; dashed 
line: Hernquist Model; dotted line: NFW Model. For each profile, the scale radius is equal to unity and the 
normalization is such that M(r = 1) = 1 
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Creation of a new bin ; P|= 2P| 

Q = Ps = Ps = 1 Ps = 2 Ps = 3 Ps 4 4 Ps *= 5 Ps = 6 Ps = 7 Ps = 8 
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Fig. 2. — Upper panel : A multi timesteps scheme is shown. Pj is computed for each bin. This figure 
shows three different bin transitions, as discussed in the text. 1: These transitions are possible because 
Pi = Pj = Q; 2: These transitions can occur only when Pj = Pj, In these cases, the particles have to 
wait until the bins are synchronized. 3: Creation of a new bin. A particle in the fundamental bin needs to 
decrease its timestep. When the new fundamental bin is created, Pj = 2Pj for every j. Lower panel: The 
red dots indicate where a force computation is necessary. 
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Fig. 3. — A 2D representation of the division into individual cubic cells. Fade grey cells represent cells 
adjacent to the grey one. / is the distance between the particle for which we want to identifiy its closest 
neighbor and the wall, and d is the distance between its closest neighbor in its own cell. 




Fig. 4. — Calculation of the distance between the considered particle, walls, and corners of the cubic cell. In 
a 2D representation, each square cell has 8 adjacent neighbors. 
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Fig. 5. — Schematic representation of a two-body system. Vectors a\ and represent the acceleration of 
particle 1 and 2. respectively. 




Fig. 6. — Schema of the reduced system 
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